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ABSTRACT 

Gravitational lensing surveys have now become large and precise enough that the in- 
terpretation of the lensing signal in current and future surveys has to take into account 
an increasing number of theoretical limitations and observational biases. Since much 
of the lensing signal is stronger in the non-linear scales, only numerical simulations 
can reproduce accurately enough the various effects one has to take into account. This 
work is the first of a series in which all gravitational lensing corrections known so far 
will be implemented in the same set of simulations using realistic mock catalogues. 
In this first paper, we present the TCS simulation suite and compare basic statistics 
such as the second and third order convergence and shear correlation functions to 
predictions for a large range of scales and redshifts. These simple tests set the range 
of validity of our simulations. We also compute the non-Gaussian covariance matrices 
of several statistical estimators, some of them are used in the Canada France Hawaii 
Telescope Lensing Survey (CFHTLenS). From the same realizations, we construct halo 
catalogues and present a series of halo properties that are required by most galaxy 
population algorithms. 

Key words: cosmology: dark matter — weak lensing — large scale structure of 
Universe — methods: systematic 



1 INTRODUCTION 

The latest measurements of the cosm ic mi- 
crowave background ( C MB) (|jarosik at al.l I2OIII : 



iPlanck Collaboration et alj I2OIII) a nd of large scale 



gala^xy surveys jYork et al 



I2OOOI : IColless et all l2003l : 

ISemboloni et all I2OO6I ) are accurate enough to constrain 
most of the cosmological parameters at the few per- 
cent level. Observations seem to converge towards a 
standard model of cosmology, in which the Universe is 
mainly filled with a uniform dark energy component, 
and about a quarter of its energy distribution consists of 



Tegmark et al 



dark matter ( Percival et al.i 2001: Eiscng t ein ct al. 2005 



20061: |Perciva!l et al.lBo07rikomatsu et alj 



20 111 : I Benjamin et al.ll2007l ). The current knowledge about 



this dark sector is so limited that an international effort has 
been set forward in order combine the best techniques and 
optimize the global c onstraints on dark energy parameters 
l|Albrecht et all l2006h . Next generation surveys including 



LSST llLSST Science Collaborations et al.ll2009l b EUCLID 
l|Beaulieu et al.ll201Gl ). SKA l|Laziol I2OO8I ), Pan-STARR^|, 
VST-KiD^I, DeH" are designed to have very high quahty 
data and large statistics, such that systematics and sec- 
ondary elTects need to be understood at the sub-percent 
level. 

Weak lensing analyses, which are based on the mea- 
surement of the degree of deformation caused by foreground 
lenses on background light sources, are particularly praised 
to detect dark matter structures. The signal allows us to 
characterize the average mass profile of foreground lenses, 
which can consist of galaxies, groups or clusters of differ- 
ent type, redshift, morphology and color, typically centered 
on a dark matter halo. The signal from the 2-point cosmic 
shear depends on seven cosmological parameters, and is es- 
pecially powerful at constraining a combination of the nor- 
malization of the matter power spectrum cts and the matter 
density Q.m. The degeneracy is further broken with mea- 
surements of the skewness and other higher-order statis- 
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tics (|Bernardeau et al]|l997l ). In the context of the global 
dark energy effort, high precision measurements of these two 
parameters are complimentary to other techniques (Bary- 
onic Acoustic Oscillations, Type lA Supernovae and cluster 
growth). It has recently been realized that weak lensing is 
also a standalone probe of the dark energy equation of state, 
through dependencies on the redshift-distanc e relation, the 
growth factor and the non-linear clusterin g (lHutereill2003 : 
lAlbrecht et~al]|2006l : iHoekstra fc JainlbOOSl ) . It is thus of the 
utmost importance to minimize the statistical, systematic 
and theoretical uncertainties associated with the signal. 

The accuracy at which one can model this signal de- 
pends on the reliability of our understanding of the lenses 
and the sources, which in turn depends on the underly- 
ing matter density field, on the specific calculation of the 
light propagation, on the exact details of the galaxy pop- 
ulation algorithm, and on our understanding of secondary 
effects. Early generatio ns of calculations were perform ed 
within linear theory (see iBartelmann fc Schneiderll200ll . for 
a review), which is known to underestimate the amount 
of structure in a large dynamical range. The projection of 
the density field on to a light cone causes the non-linear 
structures of the matter field to contribute to angles of 
a few tens of a rc minutes for sources at redshift z — 1 
II Jain et al.ll2000l ). Better r esults were obtained fro m higher- 
order p erturbation theory (iBernardeau et al.ll2002l ') , the halo 
model JCoorav fc Sheth' '20021) or with non-linear fitting 
formula jSmith et al. 2003). However these models gener- 



mology. The 'L-series' consists of 16 realizations that cover 
to lower-A; modes only, the 'H-series' covers the quasi-linear 
regime, and only one realization res olves the k 1 fcMpc ~^ 



obtained from N-body simulations (iBlandford et al 




1991 


Premadi et al.l 19981: White fc Hul I2OO0I: Jain et al. 




200c 


Hilbert et al.ll2009l: Sato et al.ll200g|: iKiessling et al.ll201ll). 



In addition, the description of the statistical uncer- 
tainty associated with weak lensing measurements needs 
to be as accurate as the signal itself, for the sake of ro- 
bustness and constraining performance. The non-linear na- 
ture of the density field on small scales tends to cor- 
relate measurements that are otherwise independent. For 
instance, the Fourier modes of the density fields grow 
independently in the linear regime, but couple together 
in the trans-linear regime (|Coles fc Chiangj |200G| ). which 
gives r ise to non-Gauss i an fe a tures in the power spec- 
trum llMeiksin fc Whiti 1 19991: iRimes fc HamiltonI l2005l: 
Nevrinck et al.l I2OO6I: iTakahashi et all l2009l : iNgan et all 
20 111 : Harnois-Deraps fc Penll2011 ). These non-Gaussian sig- 



natures propagate in weak lensing measurements, as ob- 
served in simulat ions (|Takada fc Jairj 20091: Pore et al.ll2009l : 
ISato et al.ll201ll ') and in the data (|Lee fc PenI |200^ . and 
have an effect on the power at which the lensing sign al can 
constrain cosmology (jPore et al.ll2009l : ILu et al.ll2010l ). The 
first goal of this paper is to provide robust non-Gaussian 
error bars on estimators commonly used in weak lensing 
analyses, including an accurate measurement of their non- 
linear covariance matrices. To match the resolution of mod- 
ern surveys, our measurements must be accurate at the sub- 
arcminute level, thus deeply probing the non-linear regime. 

One of the main limitation of the existing simula- 
tions is either that they are not resolving small enough 
scales, or they are limited in terms of number of realiza- 
tions. For example, t he Coyote Universe simulation suite 
l|Lawrence et al.l I2OI0I ) used cosmological volumes of 1300 

Mpc per side, organized in three series for each cos- 



scales. Also, the analyses carried bv iKiessling et al 



2011 ) 



20111 ) 



was based on the SUNGLASS pipeline ( Kiessling et al.l 
which relies on the GADGET2 tree-PM code. They used 
100 simulations with 512'^ particles and a box size of 512 
h^^Mpc. Running a new series is necessary for no existing 
suite has both the statistics and the resolution required to 
extract estimates of the non-Gaussian covariance matrix on 
2- and 3-points functions with sub arc-minute accuracy. 

In this work, we are constructing the TCS simulation 
suite with CUBEP3M with eight times more particles than 
the above mentioned SUNGLASS series, in a volume more 
than forty times smaller, thus probing the non-linear regime 
much deeper. The resources that GADGET2 would require 
for our task is very larg^fl, which explains why we opted for 
an alternative N-body code. 

The second goal of this paper is to set the stage to 
start quantifying in details the main weak lensing secondary 
effects, which are often neglected, overlooked, and in any 
case have never been included all at once in the same simu- 
lated mock catalogues. With the forecasted accuracy of the 
next generation surveys, these need to be carefully exam- 
ined, since they are likely to contribute to a large portion 
of the theoretical uncertainty. For instance, the impact of 
intrinsic alignment nee ds to be quantified in order to cali- 
brate the lensing signal l|Hevmans et al.ll2004l ). This effect is 
caused by the fact that galaxies that live in a same cluster 
are subject to a coherent tidal force, which tends to com- 
press them along the direction t o the centre of mass of th e 
halo H eavens fc Peacockl 1 19881 : ISchneider fc Bridld bOloD . 
Another secondary effect that needs to be examined is the 
so-called intrinsic alignment-lensing correlation (sometimes 
referred to shear-ellipticity correlation or contamination), a 
correlation that exists between the intrinsic alignment of the 
foreground galaxies and the shear signal of the same galax- 
ies on background sources, provided the orientation of the 
foreground galajcy correlates wit h the tidal field it is sub- 
jected to (|Hirata fc SeliakI |2004| ). Source clustering is an- 
other important secondary effect, which is caused by the 
fact that sources are not uniformly distributed : regions of 
the sky with more s ources are likely to provide a stronger 
weak lensing signal (|BernardeaiJll99 j ). Also to be tested is 
the possible intrinsic alignment of galaxies with voids, an 
effect which was previ ously found to be consistent with zero 
iHevmans et al .' 2006). 

Many secondary effects have already been studied 
iHevmans et al.ll2006l : ISemboloni et al.ll2008l ), but the sta- 
tistical accuracy and the resolution were limited. These pre- 
vious works need to be extended. More importantly, these 
have been studied separately, and we do not yet understand 
how they blend together. The only way to accurately quan- 
tify their impact in the data is by measuring their combined 
contribution in accurate mock galaxy catalogues. There ex- 
ists a vari ety of analytical methods to generate mock matter 
densities (|Coles fc Joneslll99ll : [Scoccimarro fc ShethI I2OO2I ) 



* For an estimator with N bins, the associated covariance matrix 
has A'^^ entries, and will have converged if estimated with at least 
realizations. 
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which one can populate with galaxies, however these are 
incapable of reproducing accurately enough the non-linear 
dynamics, predominant at the galactic scale. N-body sim- 
ulations, on the other hand, are perfectly suited for this 
exercise, as one can test separately the accuracy of the un- 
derlying density fields, the halofinder algorithm, the galaxy 
pop ulation scheme and the prop osed weak lensing estimator 
(see I Forero- Romero et all |2007|) for example). We thus plan 
to construct a large sample of mock galaxy catalogues with, 
in mind, to quantify the uncertainty on these secondary ef- 
fects in a later stage. 

This paper is addressing the first step in the construc- 
tion, which is the determination of key properties of the 
underlying dark matter haloes, in the cosmological context 
under study. On the longer term, our catalogues can also be 
used to test new ideas that might contribute to the system- 
atics of weak lensing signals. We do not attempt to measure 
all these effects here, but rather focus on the creation of the 
gravitational lenses and their associated halo catalogues, in 
order to set the scene for future work and know the lim- 
itations of our simulations. An interesting analysis, in the 
sa me spirit as that which we plan to perform, was carried 
in ISemboloni et al] (|201ll ). however we are improving both 
in statistics and resolution, plus we are extracting the co- 
variance matrices on the estimators, and our objective is to 
include a number of secondary effects that were overlooked. 

In this paper, however, we do not construct simula- 
tions that allows the study of baryonic feedback on the dark 
matter distribution. Recent work seems to suggest that ef- 
fects such as AGN feedback and supernovae winds could 
impact noticeably the matter distribution in the Universe 
jSemboloni et al.i l201l[ ). Although this is something that 
could be significant for the interpretation of the lensing sig- 
nal, especially at small angular scales, ignoring it does not 
diminish the value of studying the combined impact of the 
other lensing and systematics secondary effects mentioned 
before in the context of pure dark matter simulations. This 
is the path we have chosen to follow in our work. One could 
imagine that future generations of simulations could imple- 
ment all the effects we are discussing here plus the effect of 
baryonic feedback. 

This paper is organized a follow. In Section[21 we briefiy 
review the theoretical background relevant for weak lensing 
studies, then we describe the design strategy of our simu- 
lations, as well as our numerical methods in Section [31 In 
Section 21 we test the performance of our simulated lines- 
of-sight. We present the weak lensing estimators and their 
non-Gaussian uncertainty in Sections [SI and [51 In Section [51 
we describe the halo catalogues, and conclude in Section [^ 



where the matrix »I'ij(0) encapsulates the distortion of the 
two dimensional image. It is determined by four components, 
namely a convergence k, two shear components 71 and 72 
that combine together into a complex shear 7 = 71 + 172, 
plus an asymmetric factor uj. The latter comes from com- 
puting the distortion matrix to second order, but even then 
it has negligible contributi on in realistic situatio ns (see, for 
example, the Appendix of ([Schneider et al.ll 19981) ). hence we 
drop it. We thus write 



K -I- 71 
72 



72 



71 



(2) 



All of these elements are locally determined by the Newto- 
nian potential $ via; 

$,11 +$,22 $.11 - $.22 ^ ,„^ 

K= ,71 = — ^ — ,72=^,12 (3) 

where the ', i' refer to derivatives with respect to the coordi- 
nate i. For a source located at a comoving distance Xa, the 
projected distortion is computed as: 



DiXs) 



(4) 



where c is the speed of light, and the comoving angular di- 
ameter distance D{x) depends on the curvature: 

D{x) = I X for A' = (5) 

[-A(~7)sin(-A'^x) for A < 

with 



(6) 



Ho is Hubble's parameter, Qm and Qa are respectively the 
ratio of the mass and dark energy densities to the critical 
density. 

The convergence field is particularly interesting theo- 
retically since it relates, through Poisson's equation, to the 
density contrast 5: 



with S is defined as: 
5(x) = P^'^l-P 



(7) 



(8) 



where p is the average matter density in the Universe and 
p(x) is the local density. In the absence of curvature, the 
comoving angular diameter distance D{x) takes the simplest 
form, and we can extract the projected convergence n up to 
a distance Xs as: 

K{e,xs)^ r wix)S{x,e)dx (9) 



2 THEORY OF WEAK LENSING 

In this section, we briefiy review the theoretical background 
relevant for weak lensing numerical analyses. The propaga- 
tion of a photon bundle emitted from a source located at f3 
and observed at an angle G in the sky is characterized by a 
Jacobian matrix A{9): 



A,, (6) 



dp 



(1) 



where W{x) is defined as 



and 



abc) = X / dx'n{x') 
The latter reduces to x(l 



X - X 



(10) 



(11) 



when the source galaxies 
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are taken to reside in a single plane at Xs- Otherwise, the 
source distribution n(x) is taken to match observed surveys. 
In this paper, however, we use the single plane case both for 
illustrative purposes, and to disentangle cleanly the sources 
from the lenses. 

From the convergence map, finally, one can extract the 
shear fi eld by solving for the gravitational potential from 
[Eq.E] H aiser fc Sauires|[T993 n In this process, one must 
be careful about the method used to perform this calcu- 
lation, since the Universe is not periodic, while simula- 
tions usually are. The edge effects can therefore contaminate 
the calculations, hence it is necessary to somehow pad the 
boundaries. 



3 NUMERICAL METHOD 

As mentioned in the introductory section, gravity is a non- 
linear process, hence the predictions from the linear theory 
of large scale structures are only valid on the largest scales. 
In the context of the detection of weak lensing, photons tra- 
jectories are probing a broad dynamical range, and are in- 
deed sensitive to galactic scale structures, where the matter 
fields are highly non-linear. Although higher-order perturba- 
tion theory can describe such systems, the accuracy of the 
calculations are limited by the complex dynamics. We thus 
rely on N-body simulations to generate reliable non-linear 
densities, then we propagate photons and extract the re- 
sulting image distortion matrix. In this section, we describe 
some of the considerations one must keep in mind when per- 
forming such calculations. 

3.1 Tiling the lenses 

Weak lensing simulations need to be optimized as a func- 
tion of the specific measurements under study. In the current 
paper, we attempt to estimate the covariance matrix for a 
number of estimators, hence we need a large number of real- 
izations, with a sub-arcminute accuracy. In an ideal world, 
one would simulate the complete past light cone that con- 
nects the observer to the light sources, for a given opening 
angle and pixel resolution. Unfortunately, for sources that 
extend to redshift of a few, this cannot be simulated all at 
once, since the far end of the cosmological volume is at an 
earlier time than the close end. It is, however, the only way 
one can model the largest radial modes of a survey. Luckily, 
these modes c ontribute very little to the weak lensing signal 
l|Limbeil 19531 ). The coherence scales of the largest structures 
which contribute to the signal rarely extends over more than 
a few times the size of large clusters, so simulation box sizes 
of the order of a few hundreds of /i~^Mpc's generally suffice 
to model the relevant structures. These simulated boxes can 
then be stacked as to create a pencil shaped volume, or a 
line-of-sight (LOS), inside of which photons are propagated. 

One can use a different s i mulat ion for each redshift 
box, as done by (| White fc Hull2000l '). but this method is 



° For this operation, we are working under a flat sky approxi- 
mation, which allows us to perform the Fourier transforms in the 
traditional plane wave basis, and to simplify the derivatives as a 
simple flnitc difference numerical calculation. 



cpu consuming, since a single LOS involves running between 
10 and 40 N-body simulations. For covariance matrix mea- 
surements, we need hundreds of these high precisions LOS, 
hence we opted for the common work around, which con- 
sists in treating different redshift density dumps of a single 
simulation as different sub-volumes of the same past light 
cone. Because the large scale structures evolve across red- 
shift slices, there exists a systematic correlation between the 
lenses that can be minimized by randomly rotating the boxes 
and shifting the origin. This procedure allows us to reduce 
by at least an order of magnitude the number of simulations 
required for our measurements. 

The next step consists in calculating the photon 
geodesies through the large scale structures, and to com- 
pute the cumulative deformation acquired along each tra- 
jectory. The most accurate calculations are performed by 
computing these geodesies in three dimensions, along their 
trajectory, starting from the observer's camera an d progress- 
ing towards higher redshifts jVale fc Whitell2003h . Such ray 
tracing methods provide one of the most reliable measure- 
ment of the cumulative shear, convergence and deflection 
angle measured at each pixel of the observer's camera, but 
need to be calculated at run time, or otherwise require 
to store the full density contrasts in memory. Cosmologi- 
cal codes that perform ray-tracing calculations at run time 
(|Kiessling et al.lboill ') typically run much slower, especially 
in high resolution of large volume simulation s, and analyses 
that use the full three-dim ensional densities (|Vale fc White! 
l2003l : iHilbert et al.l [2OO9I ) have a large memory footprint, 
two limiting factors for the task at hand here. 

However, it was shown l|Vale fc Whitell2003l ') that differ- 
ences in lensing maps obtained from the mid-plane 'tiling' 
technique is less than 0.1 per cent, with nearly indistin- 
guishable effects on the two- and three-point functions. 
This tiling approach consists in collapsing the cosmological 
sub- volumes into their mid-planes, creating two-dimensional 
slabs - or tiles - and calculating the geodesies on these thin 
lenses. Typically, all tiles have the same comoving dimen- 
sion, and the past light cone is interpolated on to a set of pix- 
els, whose sizes correspond to the angular resolution of the 
simulated telescope. In the weak lensing regime, these trajec- 
tories are close to straight line s, such that Born's a pproxima- 
tion is very accurate ( Sch neider et al.lll998l : FVale fc Whit3 
l2003h . In this paper, we thus opt for a line-of-sight integra- 
tion along the unperturbed photon paths. We nevertheless 
stored the full simulated lenses, allow ing future analysis to 
test how m uch post-Born calculations (^ Schneider et al.lll998l : 
iKrause fc^Hirata 2010 ) affect the results. 

Computations in this setting are fast, but not optimal: 
the interpolation becomes very strong at low redshift, thus 
increasing the impact of the simulation softening length. 
More over, a large portion of the simulated volume is left 
unused: the past light cone, shaped like a truncated pyra- 
mid, is extracted from a cuboid. We could in principle shuffle 
the projections axis and the origin once more to create addi- 
tional LOS. However, this would inevitably produce a small 
amount of extra correlation between different realizations, 
that would propagate and contaminate the covariance ma- 
trix with extra non-Gaussian features. On the other hand, 
this would allow for hundred of additional, weakly corre- 
lated, LOS, which might be of use for future analyses that 
require even larger statistics. 
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3.2 N-Body simulations 



The N-body simulations are prod uced by CUBEP3M, an 
improved version of PMFAST (M erz et al.l [2005 ) that is 
both MPI and OPENMP parallel, memory local and also 
allows for particle-particle (pp) interaction at the sub- 
grid level. 1024^ particles are placed on a 2048'^ grid and 
their initial grid displacements and velociti es are calcu- 
lated from the Zel'dovich ap proximation (,ZerDovichlll970l : 
IShandarin fc ZeldovichI ligsgD with a transfer function ob- 
tained from CAMB (|Seliak fc Zaldarriagal Il996l '). The fol- 
lowing cosmological parameters are used both as simula- 
tion input and in our theoretical predictions : Qa = 0.721, 
= 0.279, ^6 = 0.046, Us = 0.96, erg = 0.817 and 
h = 0.701. 

As mentioned in the introduction, this work is meant 
to outperform the dynamical range of previous weak lensing 
simulations: we need sub-arc minute precision and a field of 
view of a few degrees per side. We designed our LOS such 
that each pixel has an opening angle of 0.21 arc minute 
on each side, with Upix = 1024^ pixels in total, for a total 
opening angle of 3.58° per side. 

In order to reduce the wasted cosmological volume that 
falls outside the past light cone, we produce d two sets of 
simu lations, somehow following the strategy of (|White fc Hul 
- which used six box sizes to z = 1. It would be compu- 
tationally too expensive to run that many distinct volumes, 
but we find that two box sizes offers a good trade off. High 
redshifts {z > 1.0) volumes are produced from simulations 
with a comoving side of L = 231.1 /i^^Mpc, while the low 
redshift ones are L = 147.0 /i~^Mpc per side. These vol- 
umes are chosen such that the boundaries of the past light 
cone intersect with the edges of the smaller box exactly at 
2: = 1 (in the given cosmology). The cone starts out at z — 0, 
then enters the larger volume, and meets its boundary at 
z = 2.0 (see Fig. [T]). Some of the outer ray bundles even- 
tually leave the simulated volume at larger redshifts larger 
than 2.0, in which case we enforce the periodicity of the sim- 
ulations. This situation applies only to the last four lenses, 
hence the total amount of repeated structures is very small. 
This is even further suppressed by the lensing kernel, which 
favors redshifts closer to z = 1 — 1.5, and by the fact such 
high redshifts have much less galaxies to start witljfl. 

With these choices of cosmological parameters and sim- 
ulation volumes, the particle's mass in the large and small 
boxes are of 1.2759 x 10^ and 3.2837 x 10* Mq respectively, 
while the comoving sub-grid softening lengths rsoft are of 
112.8 and 71.8 h'^kp^ 

The initial redshifts are selected such as to optimize 
both the run time and the accuracy of the N-body code. 
These are chosen to be Zi = 40.0 and 200.0 for the large 
and small box respectively. The reason for selecting differ- 
ent starting redshifts resides in the fact that the smaller box 



° We have used some of the otherwise wasted volume on occasions 
to extract lenses at low redshifts for some simulations that could 
not run until the end (reaching hard walltime limit, computer 
cluster downtime, etc.). In that process, we made sure that there 
were no overlap with the original lens. 

^aoft is defined to be a tenth of a grid cell, and is enforced by 
a sharp cutoff in the force of gravity for particle pairs separated 
by smaller distances. 




231.1 Mpc/h 



147 Mpc/h 



Figure 1. Geometry of the lines-of-sight. The global simulated 
volume consists of two adjacent rectangular prisms, collapsed as a 
series of thin lensed. As explained in the text, high redshift lenses 
have larger comoving volume, but the same number of grid cells, 
or pixels; this is meant to reduce the volume that falls outside 
of the past light cone. The observer sits at z = 0; the junction 
between the small (lower-z) and large (higher-z) simulation boxes 
occurs at 2 = 1; the past light cone escapes the simulated volume 
beyond z = 2, and we exploit the periodicity of the boundary 
condition to populate the edges of the most remote lenses and 
halo catalogues; we store lenses and haloes up to z = 3. 



Table 1. Redshifts of the lenses. The projections for zi > 1.0 are 
produced with L = 231.1 /i~^Mpc simulations, while those for 
lower zi are obtained from L = 147.0 /i~^Mpc, as described in 



the text. 


3.004 
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1.728 
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1.371 


1.215 


1.071 


0.961 


0.881 


0.804 


0.730 


0.659 


0.591 


0.526 


0.463 


0.402 


0.344 


0.287 


0.232 


0.178 


0.126 


0.075 


0.025 







is probing smaller scales, hence it needs to start earlier, at a 
time where the Nyquist frequency of our grid is well in the 
linear regime. Each simulation is then evolved with cubepSm 
on 8 nodes of the Tightly Cou pled System superc omputer at 
the SciNet HPC Consortium i|Loken et a"Lll2010l ') - to which 
system we dedicate the name of the simulation suite. At each 
of the redshifts tabulated in Table [U the dark matter parti- 
cles are placed on t o a 2048^ grid using a 'cloud -in-cell' in- 
terpolation scheme (|Hocknev fc Eastwood|[T98lh . That grid 
is then collapsed in three different ways - along each of the 
3 cartesian axes - producing triplets of slabs. These lens 
redshifts, zi, are found by slicing into cubes our simulated 
volume, starting and ending ai z — 0.0 and z = 3.0 re- 
spectively, and solving for the redshift at the centre of the 
comoving box. 

With this configuration, we solve [Eq.[U] numerically for 
each pixel. We convert the x integral into a discrete sum at 
the lens locations xi^i)- The infinitesimal element dx be- 
comes dL/ugrid, where Ug^id = 2048 and L = 147.0 or 231.1 
/i~^Mpc, depending on the redshift of the lens. Under the 
single source plane approximation, we can thus write the 
convergence field as 



k(x) = 



2c2 



J2 + ^Ox(zO(l - X{zi)/X{^s))dx (12) 



where 5(x) is the two-dimensional density contrast. 

To avoid edge effects when computing the shear fields. 
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we perform the Fourier transforms on the full periodic slabs, 
i.e. before the interpolation on to the lenses. 



4 TESTING THE SIMULATIONS 

In this section, we quantify the resolution and accuracy of 
the simulated weak lensing simulations. We first measure 
the power spectrum of the simulated three dimensional den- 
sity fields, i.e. before the collapse and pixel interpolation, 
and compar e to t he non-linear theoretical predictions of 
ISmith et al.l (|2003l i^. We then estimate the angular power 
spectrum of the simulated lines-of-sight, compare to non- 
linear predictions, and extract the effective resolution of the 
simulated fields. 



4.1 Dark Matter Power Spectrum 

The power spectrum of matter density P{k) is a fast and 
informative test of the quality of the simulations. It probes 
the growth of structures at all scale available within the sim- 
ulations, and comparison with a reliable theoretical models 
informs us of the accuracy and the resolution limit. For a 
given density contrast 5(x), the power spectrum can be cal- 
culated from its Fourier transform 5(k) as: 

(5(k)5(k')) = (27r)'fo(k - k')P(k) (13) 

where the angle brackets refer a volume average, and the 
Dirac delta function selects identical Fourier modes. We 
present the power spectrum for our 185 simulations for two 
redshifts, z = 0.961 and z — 0.025, in Fig. [2] The error bars 
are the la deviation in the sampling. 

Resolution limits are mainly determined by the soften- 
ing length in the gravitational force, and cause an abrupt 
drop in the observed power spectrum in small scales. In our 
simulations, we observe that the structures seem to be well 
model down to fc = 20.0/iMpc~'^, which corresponds to a 
comoving length of about 315 kpc. This d rop of power 
can be modeled, following (|Vale fc White|[20o3 ) . by a Gaus- 
sian filtering of the form exp[— fc^cTg] in the power spectrum, 
where cTg — 0.155L/Ngrid- This filter is in comoving coor- 
dinates, hence it is straight forward to calculate the angle 
subtended by lag at z 1.0, the redshift that dominates 
the lensing kernel. From this rough estimate, we find that 
the softening angle is Osoft = 0.148 arcsec. We note that this 
technique depends on the details of the N-body code, and 
that a more accurate estimate of the softening length can 
be found from a comparison of the angular power spectrum 
with non- linear theory, which we post-pone until section 
Otherwise, we also observe that the higher redshift power 
is slighly on the low-side of predictions for 0.4 < k < 1.0 
hMpc'^ , while the lower redshift power is sligthly higher 
for k > 1.0 hMpc-^. These can be caused by a number of 
effects, from finite box size to residual uncertainty from nu- 
merical integration, and inevitably propagate in the calcula- 
tions of past light cone. However, deviations from HALOFIT 



* The version of HALOFIT that was used in our calculations does 
not incorporate recent corrections that were made in December 
2011 




10"' 10° io' 



k[hlMpc\ 

Figure 2. Power spectrum of 185 N-body simulations, at red- 
shifts of 0.961 (bottom curve) and 0.025 (top curve). The solid 
and dashed lines are the non-linear predictions, with and without 
the Gaussian filter. The error bars shown here are the standard 
deviation over our sampling. We observe a slight over estimate of 
power in the simulations for scales smaller than k = 3.0/iMpc~^. 
This is caused by a known loss of accuracy in the predictions, as 
one progresses deep in the non-linear regime. 

are observed in most N-body codes, and do not affect qual- 
itatively the final results, as long as internal consistency is 
preserved. 

In linear theory of structure formation, different Fourier 
modes of the matter density grow independently, such that 
the error bars on the power spectrum are well described 
by Gaussian statistics. For non-linear scales however, the 
phases of different Fourier modes start to couple together 
jMeiksin fc Whitell999l : IColes fc Chian3l2000l : IChiang et all 
'2002*), hence the two-point function no longer contains all 
the information about the fields. Higher order statistics, i.e. 
bispectrum, trispectrum, are then needed in order to com- 
pletely characterize the field. Theoretical calculations can 
describe these quantities with a reasonable accuracy in the 
translinear regime, but N-body simulations provide the most 
accurate estimates at all scales, as long as convergence is 
achieved. For a covariance matrix with A^^ bins, one needs 
at least the same amount of simulations. Our plan is to or- 
ganize the lensing estimators in about ten bins, hence 185 
simulations are enough. We recall here that the simulations 
described in this paper are primarily used to estimate the 
non-Gaussian covariance matrices of various weak lensing 
estimators, hence it is essential to pin down the source non- 
Gaussian features, and to monitor how these propagate in 
our calculations. 

The power spectrum covariance matrix is defined as 

C(fc,fc') = (P(fc) - P{k)){P{k') - P{k')) (14) 

where the over-bar refers to the best estimate of the mean. 
The amount of correlation between different scales is bet- 
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Figure 3. Normalized covariance matrix of tiie density power 
spectrum, measured from of 185 N-body simulations, at redshift 
of 0.961. Modes at fc ~ 0.5/iMpc~^, corresponding to 9 ~ 18.35 
arcmin, are more than 40 per cent correlated. 



ter visualized with the cross-correlation coefficient matrix, 
which is obtained from C{k,k') via 



p{k,k') = 



C{k,k') 



(15) 



^C{k,k)C{k',k') 

and is shown for z — 0.961 in Fig. [3] We see that it is al- 
most diagonal in larger scales (lower k), while measurements 
become correlated as we progress towards smaller scales 
(higher fc). At fc 0.5/iMpc~^, for instance, the Fourier 
modes are intrinsically more that 40 per cent correlated. 
This corresponds to an angle of ~ 18.35 arcmin on the 
sky, and I ~ 1180. 

We stress that this correlation is intrinsic to the den- 
sity fields, and that extracting this covariance matrix from 
an actual galaxy survey is a delicate task. Recent results 
have shown that neglecting either the non-Gaussian nature 
of the uncertainty, or the survey selection function can sig- 
nificantly underestimate the error on the powe r spectrum, 
even at scales traditionally conside red as linear (jNgan et al.l 
I2OIII : iHarnois-Deraps fc Penlboill ) . 

4.2 Angular Power 

In order to quantify the resolution of our lensing maps, we 
measure the angular power spectrum of the k{0) fields, and 
compare the results with the non-linear predictions from 
ISmith et all ([2003'). The power spectrum of the convergence 
field is defined as: 



(k(^i)k(£2)) = (27r)'fc(£i -f £2)P«(^i) 



(16) 



where £ is the Fourier component corresponding to the real 
space vector G. and where, again, the angle brackets refer 
to an angle average. The convergence power spectrum, esti- 




Figure 4. Convergence power spectrum, measured from 185 N- 
body simulations, where the source redshift distribution is a Dirac 
delta function at 2 ~ 3.0, 1.5, 1.0 and 0.5 (top to bottom sym- 
bols). The solid lines corre spond to the non-hnear predictions 
llSeliak fc Zaldarriaga|[l996l) . the linear prediction at = 3 is 
represented by the dashed hne, and the error bars are the la 
standard deviation over our sampling. We observe a slight over- 
estimate of power in the simulations for z = 3 and £ > 1000 com- 
pared to non-linear predictions (sohd hne), and a more important 
bias for lower redshifts. This roots in the lack of predicted power 
which is also visible in the smallest scales of the three dimensional 
dark matter power spectrum (i.e. Fig. [2]l. Similar trends are ob- 
served in the Coyote Universe and SUNGLASS simulation suites. 
The low-£ power seems also to be in excess in the simulations, 
however predictions are still with the error bars. 



mated from our simulations, is shown in Fig.|4]where the er- 
ror bars are the la standard deviation on the sampling. It is 
presented in the dimensionless form, i.e. ^(£-1- 1) / (2tv)Pk (^ 
which has the advantage of showing clearly which angles (or 
multipoles) are probing the linear regime. When compared 
with the non-linear theoretical model, we find a good agree- 
ment in the linear regime, while the theoretical predictions 
slightly underestimate th e power for £ > 100 0, consistent 
with the observations of iHilbert et al.l l|2009l ). The strong 
drop at ^ ~ 30000 is caused by limitations in the resolution, 
which corresponds to an angle of about 0.7 arcmin. 

As mentioned earlier, the smallest angles of weak lens- 
ing observations are probing the non-linear regime of the 
underlying density field, and it is known that the statistics 
describing the uncertaint y on the weak len sing power spec- 
trum are non-Gaussian ( Pore et al.ll2009l ). Although most 
of the departures form Gaussianity in the data are currently 
lost in the observation noise, future lensing surveys are ex- 
pected to improve enough on statistics and systematics such 
that non-Gaussian features will become significant. The non- 



^ In this paper, wc use two notations interchangeably to denote 
the convergence angular power spectrum : Pk(<?) and Ci. 
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Figure 5. Cross-correlation coefficient matrix of the (dimension- 
full) convergence power spectrum, measured from 185 LOS. We 
observe a strong correlation for £ of a few thousand. 

linear dynamics effectively correlate the error bars in small 
scales, as seen in Fig. [S] As expected, we observe that all 
the multipoles with £ > 1000 are more than 40 per cent 
correlated. 



dom points and construct mock catalogues, from which we 
extract the 2-point correlation function measurements. The 
positions are purely random within the 12.84 deg^ patches, 
and the values at each point are interpolated from the simu- 
lated grids. This completely bypasses the galaxy population 
algorithm, which will be addressed in future work. 

5.1 Shear 

In current cosmic shear analyses, one of the strongest weak 
lensing signal comes from a measurement of the 2-point cor- 
relation function in the shear of galaxy, which is defined as: 

C,,(e) = (7,(0')7.(e + 0')> (17) 

where i and j refer to a pair of galaxies separated by angle 

6 — In the absence of gravitational lenses, (,ij{9) averages 
out to zero, hence a positive signal indicates a detection of 
cosmic shear. 

Since the intrinsic distortion produced by a single mas- 
sive object is exclusively aligned in the tangential direction 
around the center of mass, it is possible to improve the lens- 
ing measurement by considering a local coordinate system, 
in which the tangential and rotated axes are defined as the 
direction perpendicular and parallel to the line joining two 
galaxies. The new, local, complex shear can be written as 

7 = 7t + i'yr, where t and r refer to tangential and radial 
shear respectively. The corresponding correlation function 
^tt and ^r-r are defined as the weighted average of the tan- 
gential and rotated shears for pairs of galaxies separated by 
an angle 6=\xi — Xj\, namely: 

_ T.WiWj'yt{xi)"it{xj) 



5 2-POINT CORRELATION FUNCTIONS 

Detection of a robust weak lensing signal from the data is a 
challenging task in itself, for many secondary effects need to 
be filtered out. The 2- and 3-point functions of the observed 
matter field are known to provide a wealth of information 
about the lensing signal. Generally, different statistical es- 
timators and filtering techniques are sensitive to different 
scales, systematics and secondary effects, and correlate the 

measurements in a unique way^ 

It was recently shown (|Vafaei et al.l 120101 ) that the 
optimal approach for measurements involving the cosmic 
shear and convergence depends on the survey geometry and 
the cosmological parameters investigated. For instance, the 
shear 2-point correlation function minimizes the correla- 
tion across different angles, while the mass aperture window 
statistics correlate different angular scales even in the ab- 
sence of signal. However, the latter has the advantage of be- 
ing more sensitive to smaller scal es, hence it is better su ited 
for surveys of limited coverage (|Schneider et al.| Il998l) . In 
the following two sections, we give a short overview of the 
most commonly used statistical estimators and present their 
signal as detected in the simulations, compared to non-linear 
predictions. 

As mentioned in section[2l the N-body simulations serve 
to construct the convergence and shear maps from the den- 
sity fields. To mimic the actual detection from a galaxy sur- 
vey, we Poisson sample each of the maps with 100000 ran- 



T,WiWj^r{xi)'yr{xj 



(19) 



where the Wi refers to weights of objects which quan- 
tifies how well the corresponding shear is measured. The 
combination of tangential and rotated shear leads to ^±(^): 

i±{e)=£,U±irr (20) 

which are conveniently related to the convergence power 
spectrum: 

dl 



2-K 



2^ 



(21) 



(22) 



where J„{x) is the nth order first kind Bessel function. Hence 
measurements of ^rr,tt give a direct handle on many cosmo- 
logical parameters that depend on Pk.{£), (seven in total, but 
mainly erg and fim). 

We show in Fig. |6] and [3 the 2-point correlation func- 
tions ^tt and r respectively, as a function of the separation 
angle and for 4 different redshifts. The error bars are the la 
standard deviation as obtained from our 185 lines of sight. 
The agreements between the simulations and the theoret- 
ical predictions are well within the error bars down to 0.6 
arcmin, which allows us to conclude that the signals are well 
resolved at least in that range. 

We next show in Fig. [8] the cross-correlation coefficient 
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Figure 6. Shear correlation function ^tt, computed from Poisson sampling the simulated shear maps in the local (tt, rr) coordinates, 
described in the text. The solid line shows non-linear theoretical predictions on the mean of the estimators. 
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Figure 7. ^rr component. 
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matrices related to the tt measurements, for source redshifts 
of 3.0 and 1.0. These show that the error bars across different 
angles are at least 50 per cent correlated for the highest 
redshift, and up to 80 per cent for lower redshift sources. 
This correlation becomes even stronger as the two angles 
become similar in size. The rr matrices are qualitatively 
similar, hence we do not show them here. 



5.2 Convergence 

Convergence - or magnification - signal has been suc- 
cessfully detected in recent a nalyses (|Scranton et al.l [20051 : 
iHildebrandt et al.|[2009l , I2OIII '). Although more challenging 
to detect in the data, it is worth the effort, first as it serves as 
a cross-check of the shear results, and second as it is theoret- 
ically cleaner: no (non-local) Fourier transforms are needed 
in the process of map construction. Following the procedure 
developed for the shear fields, we calculate the two-point 
K — K correlation function from the Poisson sampling of the 
simulated convergence maps. In Fig. [9l we present our re- 
sults, compared with the non-linear theoretical predictions, 
and find that the agreement extends deep under the arc 
minute at all redshifts as well. 

We have also measured the cross-correlation coefficient 
matrices corresponding to these measurements, and because 
they are very similar to the "/tt matrices, we do not present 
the matrices here, but make them available upon request. 



6 WINDOW-INTEGRATED CORRELATION 
FUNCTIONS 

As briefly explained earlier, depending on the survey, the 
2-point correlation functions are not always the best way 
to measure the cosmic shear or convergence. Window- 
integrated correlation functions, such as the mass aperture 
varian ce, give a second hand le on many cosmological param- 
eters l|Schneider et al.lll998l ). and are also used in galaxy- 
galaxy lensing and cluster lensing. Generally, the method 
consists in measuring the mean, the variance or even higher- 
order statistics of a given lensing field, which was before- 
hand convolved with either a 'top hat' or a compensated 
filter of variable smoothening angle 8. The integrated re- 
sults are computed as a function of the smoothening angle, 
and can usually be compared with theoretical predictions 
for the lensing power spectrum. 

When extracting this kind of estimator on the shear sig- 
nal, the choice of filter matters. One of the advantage of the 
top hat filter is that it probes scales as large as the field of 
view, whereas compensated filters are limited by a damping 
tail in the window shape that prevents using the boundaries 
of the patch. Another advantage of the top hat filter is that 
it yields a signal to noise tha t is optimal for a skewness mea- 
surement (| Vafaei et al.ll2010l ). On the other hand, a compen- 
sated mass aperture filter is more sensitive to smaller scales, 
hence it does a better jobs at recovering signal on surveys 
with smaller area. In add ition, it is measured directly from 
the t angential shear field l|Kaiser fc Sauires|[l993 : ISchneideij 
1 19961 ). hence is not affected by the systematic and statistical 
uncertainties involved in the reconstruction of the conver- 
gence field. For these reasons, different analyses will prefer 



different estimator, hence we consider both kinds of filters 
for completeness. 

The top hat filter is a circular aperture of radius 9, out- 
side of which the signal is cut to zero. The filtering process 
effectively measures the total shear in a filtered region of a 
map, 7, for a given opening angle 9. We repeat such mea- 
surement over many patches on each map and measure the 
top hat variance (|7|^)th, which is related to the conver- 
gence power spectrum Pk as: 

{\f{9)\}TH = ^ J dUP4£)WTHm (23) 

with Wth{19) = *fjyP (|Kaiseij|l99j ). This convolution 
process is done with Fourier transforms, and each of the 
maps are zero-padded in order to reduce the effect of bound- 
aries. We compare our measurements from the simulations 
with non-linear predictions in Fig. 1101 as a function of the 
opening angle of the top hat filter, and find a good agreement 
at all redshift, although 2 — 0.075 exhibits a small positive 
bias. The cross-correlation matrices are presented in Fig. 1111 
and show that there is a strong correlation between most 
measurements. 

We next consider a compensated aperture filter, which 
is constructed from the local tangential shear in mock cat- 
alogues, where the role of the central galaxy in the pair 
counting is played by the c entre of the filter. Th e aperture 
mass Map is then given by (|Schneider et al.lll998l l: 

Map{9) = J cFdQe{d)'^t{d) (24) 

where Q is a weight function with support £ ^,9] and 
which takes the shape: 

0«W = J.(^)(l-J) (25) 

We then calculate the variance {Map) across the map, for 
all available angles, which is also related to the convergence 
power spectrum P„ by: 

{Mi^m = — / dup,{t)Wapm (26) 

27r Jo 

where Wap{l9) = '^^^0^- We present in Fig. [H our mea- 
surements of {Map) from the simulations, as a function of 
smoothing scale 9. Over plotted are the theoretical predic- 
tions obtained from [Eq. [26] with a non-linear convergence 
power spectrum. We observe that for redshifts larger than 
one, the agreement extends down to the arc minute, whereas 
lower redshifts suffer from a lack of variance at angles of a 
few arc minutes. This is caused by limitations in the res- 
olution due to strong zooming from the simulation grid 
on to the pixel map. We recall that an aperture mass is 
constructed with a compensated filter, which has a strong 
damping tail, hence for an opening angle 9, it is really sensi- 
tive to ~ S/5. This thus approaches simulation resolution at 
very low redshifts. This drop is also expected from the top 
hat variance, but appears at much smaller smoothing an- 
gles. The cross-correlation coefficient matrices are presented 
in Fig. 1131 and show that most measurements are close to 
60 per cent correlated. The smallest angles probe scales that 
are approaching the pixel resolution, hence there is very lit- 
tle cross-correlation. 
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Figure 8. Cross-correlation coefficient matrix of the tt two-point function, with the source plane at 2 ~ 3.0 (left) and z ~ 1.0 (right). 
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Figure 9. Convergence correlation function ^kk, at four different redshifts, constructed from Poisson samplings of the simulated maps. 
The error bars are the Icr standard deviation in the 185 realizations. 



7 WINDOWED STATISTICS ON 
CONVERGENCE MAPS 

Window statistics performed directly on the convergence 
fields serve as an important test of the accuracy and pre- 
cision of the simulations, since the calculations here can be 
done directly on the grid, i.e. without the Poisson sampling. 
We smooth the K-maps with filters identical to those used 
in the last section and proceed to calculate the top hat vari- 
ance {k^{9))th and mass aperture variance {Map{9)) as a 



function of the filter opening angle. In the latter case, the 
choice of compensated filter ( i.e. the equivalent of Qb{'&) in 
[Eq.[25] ) is given by Uei'd) (ISchneider et al.lll998l ). where: 



9 



612 



(27) 
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Figure 10. Top hat variance, (|7p}, measured from shear maps, for 4 different redshifts. 
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Figure 11. Cross-correlation coefficient matrix of the top hat variance, with the source plane fixed at z ~ 3.0 (left) and z ~ 1.0 (right). 



In analogy to [Eq. [21] , the Map estimator that is obtained 
from the convergence maps is calculated as: 



(28) 



We present in Fig.[T3]and[T5]the comparison with non-linear 
predictions, and observe a signal that is almost identical 
with the corresponding shear estimator. Namely, the agree- 
ment is good at all redshifts for the top hat variance, with 
only a slight bias at the lowest redshifts, whereas the mass 
aperture variance shows a lack of signal at angles of a few 



arc mirmtes for low redshifts. Once again, this comes as no 
surprise since this approaches the grid size. 



We also show the three-point function {M'ap{9)) and 
{k^{9))th in Fig. ll6l and ll7l which allow to break the degen- 
eracy between as and ilm- We observe that our simulations 
tend to overestimate the mass aperture predictions by la at 
low angles. This comes again from the fact that the aperture 
filter is sensitive to about one fifth of the total opening angle 
probed. Hence the discrepancy observed at ~ 2 arcmin is 
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Figure 12. Aperture mass variance (M^p) measured from shear maps. The apparent discrepancy between simulations and theoretical 
predictions at low redshift is caused by resolution limits, where the small angles actually probe scales that are very close to the grid size. 
The intrinsic pixel size of 0.21 arcmin has effects up to scales of about 1.0 arcmin. 
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Figure 13. Cross-correlation coefficient matrix of the aperture window integrated shear, with the source plane at 2: ~ 3.0 (left) and 
2 ~ 1.0 (right). The correlation seems to vanish in the first bin, however, for reasons discussed in the text, these bins are under-resolved, 
hence the largest modes . 

mainly probing scales of 0.4 arcmin, "which is of the order of 8 HALO CATALOGUES 
the pixel size. 

This section briefly describes how the halo catalogues 
are created, and presents a fe"w of the properties that 
are extracted on individual objects. We recall that our 
is to construct mock galaxy catalogues on which many 
secondary effects will be quantified. As mentioned be- 
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Figure 14. Aperture mass variance (A/^p), measured directly from the convergence maps. We recall that the effect of finite pixel size is 
felt to larger angular scales - up to about one arc minute - with this estimator. 
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Figure 15. (|7P)th measured directly from the convergence maps. We see that the agreement is excellent at all redshifts. 
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Figure 16. (M^p) measured directly from the convergence maps. We recall that the effect of finite pixel size is felt to larger angular 
scales - up to about one arc minute in lower redshift lenses - with this estimator, which explains the damped tail at small angles in the 
z = 0.075 plot. 
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Figure 17. {\^\^)th measured directly from the convergence maps. 
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fore, we do not attempt to populate the haloes in this 
paper, since this is a challenge on its own, and we 
wish to factor out this problem for now. In the future, 
though, one could follow (jHevmans et al.l [20061 '). and pop- 
ulate the haloes with galaxies under the conditio nal lu- 
minosity function of (ICoorav fc Milosavlievid '2005) , then 
assign ellipticity fol l owing the elliptical o r spiral model 
l| Heavens et al.1 12OO0I : iHevmans et al. 20o3'). Al ternatively, 
one c ould use the GALI CS (|Hatton et all 120031 ') and MO- 
MAF iBlaizot et aLlbOOST ) to create mock galaxy catalogues 
directly our halo catalogues, foll owing the prescription de- 
scribed in (jForero- Romero et al.i ,2007). Our plan is to in- 
corporate, for the first time and in a systematic way, 
the prescriptions for intrinsic ali gnment, s ource clustering, 
ellipticity-shear correlati on, etc. IjPe Lucia & Blaizot 20071: 
ISchneider fc BridldbOld ) all at once. This is crucial in or- 
der to interpret correctly the signal from the data, which 
contains all these contributions. 

In this work, haloes are constructed from the mat- 
ter density fields w ith a spherical overdensity algorithm 
l|Cole fc LacevI Il996l ) . The first step is to assign the dark 
matter particles on a 2048^ grid with the Nearest Grid Point 
scheme (Hoc kney fc Ea stwood 198ll ') , and identify local den- 
sity maxima. The halo finder then ranks these halo candi- 
dates in decreasing order of peak height, and for those which 
are above an inspection threshold value, it grows a spheri- 
cal volume centered on the peak, computing for each shell 
the integrated overdensity until it drops under the predicted 
critical value. The haloes that are analyzed first are then re- 
moved from the density field, in preparation for the inspec- 
tion of lower mass candidates. This prevents particles from 
contributing multiple times, but at the same time limits the 
resolution on sub-structures of the largest haloes. This is 
a mild cost for the purpose of these catalogues, which are 
populated with low multiplicity of galaxies, and thus depend 
rather weakly on the sub-halo structures. Finally, for each 
halo, we measure the mass, the centre-of-mass (CM) and 
peak positions, the CM velocity, the velocity dispersion, the 
angular momentum in the CM frame, and the inertia matrix 
o"ij , which allows us to assign galaxy population beyond the 
sole halo mass. Although the latter quantity is biased by the 
fact that we are only searching for spherical regions, we still 
recover enough information about the shape and orientation. 

In all the plots presented in this section, we present 
properties of the haloes that populate the full simulation 
box, even though, in the final mock catalogues, we keep only 
those that sit inside the past light cone. We apply the coor- 
dinate rotation and the random shifting of the origin that 
was performed on the lensing slabs, such that the halo cat- 
alogues and the lenses trace the same underlying density 
field. 

To quantify the accuracy of the halo catalogues, we 
first extract the power spectrum of the distribution, and 
compare the results with the measurements from dark 
matter particles by computing the halo bias, defined as 
b{k) = Phaio{k) / Pparticie{k) . From Fig. [TSl we observe 
that both the shape and redsh i ft de pendence agree well 
with the results from llhev et al.l (|201ll ). We observe in Fig. 
15] that t he halo mass functi on is in good agreement with 




k\hlMpc] 



Figure 18. Halo bias, for redshifts 0.025 (bottom curve) and 
0.961 (top curve). 
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Press fc Schechter (1974) and lSheth fc Tormen (20o3) in the 



range 1 x 10 — 2 x 10 M©. Higher redshift seems to better 



Figure 19. top: Halo mass function, compared to predictions, 
for redsliift 0.025, 1.071 and 3.004 (top to bottom curves), bot- 
tom: Ration of the mass function to the theoretical predictions of 
Sheth- Tormen. 



fit Press-Schechter, while lower redshifts are better described 
by the Sheth- Tormen model. 

We now present some of the relationships between the 
main kinematic quantities of these haloes. In Fig. 1201 [2l1 and 
1221 we show, for z = 0.025, 0.961 and 3.004 respectively, and 
from top left to bottom right: the angular momentum L ver- 
sus the mass m, the halo CM velocity Vcm versus m, L versus 
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Vcm, then the dispersion velocity Vdisp versus m, L versus 
Vdisp, and finally a versus m, where a is the angle average of 
the standard deviation on the centre of mass position. It is 
related to the inertia matrix via = ^ '^ij- The mass 

is in Mq, the velocity is in km s~^, the distances in /i~^Mpc, 
and the angular momentum in /i~^Mpc km s~^. The plots 
are bivariate histograms, arranged in 50x50 bins, and the 
color bars are in logarithm scale for higher contrasts. 



9 CONCLUSION 

This paper has two principal objectives: 1) measure accurate 
non-Gaussian covariance matrices at the sub-arcminute pre- 
cision on the principal weak lensing estimators, and 2) set 
the stage for systematic studies of secondary effects, and es- 
pecially how their combination impacts the lensing signal. 
We have generated a set of 185 high resolution N-body sim- 
ulations from which we constructed past light cones with 
a ray tracing algorithm. The weak lensing signal is accu- 
rately resolved from a few degrees down to below the arc 
minute. Thanks to the large statistics, we have measured 
non- Gaussian error bars on a variety of weak lensing esti- 
mators, including 2-point correlation functions on shear and 
convergence, and window integrated estimators such as the 
mass aperture. In each case, we compared our results with 
non-linear theoretical predictions at a few redshifts and ob- 
tained a good agreement, which testifies the quality of the 
simulations. The extracted error bars are essential for a cor- 
rect estimate of many lensing signals, including sensitivity 
to cosmological parameters, which so far relied on Gaussian 
assumptions, or on numerical estimates that were not resolv- 
ing the complete dynamical range. With the next generation 
of lensing survey, non-Gaussian error bars are expected to 
deviate significantly from Gaussian prescriptions, therefore 
techniques such as those presented here will be required for 
robust estimates. 

In addition, we measured the covariance matrices of 
these weak lensing estimators, and we show that the error 
bars between different angular measurements are at least 
50 per cent correlated, with regions up to 90 per cent cor- 
related, as the two angles become closer. This needs to be 
taken into account when propagating the non-Gaussian un- 
certainty from these estimates on to other parameters. 

We also generated a series of halo mock catalogues that 
are coupled to the gravitational lenses constructed, for fu- 
ture independent studies of secondary signals and alternate 
testing of lensing estimators. Within the CFHTLenS collab- 
oration, these catalogues will be part of the CLONE project. 
We finally present an overview some key kinematic variables 
in our halo catalogues, including the mass dependence of the 
velocity dispersion and of the angular momentum. Our near 
term goal is to include effects such as intrinsic alignment, 
source clustering, etc. in our galaxy population algorithm 
and quantify their combined contribution. In addition, we 
plan to quantify the impact of post-Born calculations on the 
non-Gaussian uncertainty and on the contamination by sec- 
ondary signals. Understanding the impact of all these effects 
is essential, for such systematic bias are likely to contribute 
significantly to future surveys such as KiDS and EUCLID. 

Aspects not included in our simul ation settings are 
feedback effects (|Semboloni et all l201ll ) and dependence 



of t he covariance ma trices on the cosmological parame- 
ters (|Eifler et al.ll2009l ). While the latter can be simply ad- 
dressed by running additional simulations, the former neces- 
sitates the construction of hydrodnamical simulations that 
implement simultaneously matter clustering at large angular 
scales and a proper modelisation of feedback effects. 
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